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We present a universal approach to the investigation of the dynamics in generalized models. In 
these models the processes that are taken into account are not restricted to specific functional 
forms. Therefore a single generalized models can describe a class of systems which share a similar 
structure. Despite this generality, the proposed approach allows us to study the dynamical properties 
of generalized models efficiently in the framework of local bifurcation theory. The approach is based 
on a normalization procedure that is used to identify natural parameters of the system. The Jacobian 
in a steady state is then derived as a function of these parameters. The analytical computation of 
local bifurcations using computer algebra reveals conditions for the local asymptotic stability of 
steady states and provides certain insights on the global dynamics of the system. The proposed 
approach yields a close connection between modelling and nonlinear dynamics. We illustrate the 
investigation of generalized models by considering examples from three different disciplines of science: 
a socio-economic model of dynastic cycles in china, a model for a coupled laser system and a general 
ecological food web. 
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INTRODUCTION 



Dynamical systems are used to study phenomena from 
diverse disciplines of science such as laser physics, pop- 
ulation ecology, socio-economic studies and many more. 
The corresponding mathematical models have often the 
form of balance equations, in which the time evolution of 
the state variables is determined by gain and loss terms. 
Depending on the processes that are taken into account, 
a single state variable can be effected by several gains 
and losses. In the modeling process the modeler usu- 
ally decides first which processes are important and need 
to be included in the model. These processes determine 
the structure of the model. In the second step each of 
the terms is described by a specific mathematical func- 
tion, which can be based on theoretical reasoning or em- 
pirical evidence. In this way a specific model for the 
phenomenon under consideration is constructed. The 
investigation of specific models is a powerful approach, 
that has, in many cases, revealed interesting insights. 
However, the specific mathematical functions on which 
these models are based depend critically on the modelers 
knowledge of the system. While the modeler may have 
much information about certain processes, others may 
be known or believed to exist which are very difficult to 
quantify. 

Specific models are therefore often based on a large 
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number of assumptions. This can be illustrated very 
well by considering an example from ecology: The so- 
called Holling type-II function is regularly used to de- 
scribe the predator-prey interaction in food webs and 
food chains[lj. This function takes major biological ef- 
fects into account. But, it can not possibly capture all 
the subtle biological details that exist in nature. Such 
details may involve the formation of swarms to confuse 
predators, adaptation of the prey to high predator densi- 
ties, density dependent changes in the predator's foraging 
strategy to name a few. In models such details are of- 
ten omitted on purpose since the modeler is interested in 
general insights that do not depend on specific proper- 
ties of the system under consideration. But on the other 
hand, certain details may have a strong impact on the 
qualitative behavior of the system. For instance, it has 
been shown that even minor variations in the functional 
form of the predator-prey interaction can have a strong 
impact on the system's stability |2|, Ll ■ 

In contrast to specific functional forms, the basic struc- 
ture of the system is much easier to determine. For in- 
stance in our ecological example it is easier to say that the 
predation depends on the population density of predator 
and prey than to derive the exact functional form that 
quantifies this dependence. It can therefore be useful to 
consider generalized models which describe the structure 
of a system in terms of gains and losses but do not restrict 
these processes to specific functional forms. If a system 
is considered in which some processes are known with a 
high degree of certainty while others remain uncertain 
it can be useful to study partially generalized models in 
which some processes are described by specific functional 
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forms, while others remain general. In other cases fully 
generalized models in which all processes are modeled by 
generalized functions can be advantageous. Generalized 
models describe systems with a higher generality than 
specific models. They depend on less assumptions and 
enable the researcher to consider the system from an ab- 
stract point of view. 

Despite the generality of generalized models, it is pos- 
sible to compute the stability of a steady state in these 
models analytically. This reveals the exact relationship 
between the qualitative features of functional forms and 
the local bifurcations of steady states • In this way the 
application of generalized models enables us to investi- 
gate the local stability of systems with a high degree of 
generality. Moreover, we can draw certain conclusions on 
the global dynamics of the system under consideration. 
For instance, the presence of chaotic or quasiperiodic dy- 
namics as well as the existence of homoclinic bifurcations 
can be deduced form certain local bifurcations of higher 
codimension. In this way the presence of chaotic dynam- 
ics can be shown for a whole class of model systems 

In this paper we present a universal approach to the in- 
vestigation of generalized models based on balance equa- 
tions. We illustrate this approach using three models 
from three different fields of science: socio-economics, 
laser physics and population dynamics. These examples 
demonstrate that the proposed approach can be used to 
study a wide range of systems. 

The paper is organised as follows. We introduce our 
approach by considering a simple socio-economic model 
in Sec.^ The underlying principles of the normalization 
procedure are discussed in more detail while we study a 
coupled laser system in Sec. lIIII Finally, we show that the 
proposed approach is applicable to more involved prob- 
lems by investigating a model for general food webs in 
Sec. El 



II. A GENERALIZED MODEL OF DYNASTIC 
CYCLES 

Chinese history is characterised by repeated transitions 
from despotism to anarchy and back jfj. This periodic 
behavior is known as the dynastic cycle. As a young dy- 
nasty rises from anarchy a period of order and prosperity 
begins. As a result the population grows. At some point 
the aging dynasty can not control the growing popula- 
tion any longer. The results are poverty, crime and civil 
war which lead back to a period of anarchy. By contrast, 
European dynasties exhibit stable behavior for extended 
periods of time. 

It has been shown in 5] that the existence of dynastic 
cycles as well as steady state behavior (e.g. stable dynas- 
ties) can be predicted by a simple three-variable model 
that describes the evolution of a society of farmers X, 
bandits Y and rulers Z. The dynamics of this model have 
been investigated in [(| . In these investigations a specific 
model was used in which the interactions between the dif- 



ferent social groups were described by specific functions. 
In the following we study a generalized version of this 
model in which the interactions are described by general 
functions. 



A. Formulation of a generalized model of dynastic 
cycles 

In the model the growth of the population of farmers 
is described by a function S(X). This function repre- 
sents agricultural production, but also losses because of 
natural mortality including diseases. The population of 
farmers is reduced by crime C(X, Y) and taxes T(X, Z). 
Since the bandits profit from the crime their number is 
assumed to increase proportional to C(X, Y). In addition 
crime increases the willingness of the farmers to support 
rulers. Therefore the number of rulers also grows pro- 
portional to C(X, Y). The rulers fight crime by reducing 
the number of bandits. This effect of the rulers on the 
bandit population is described by a function L(Y 7 Z). Fi- 
nally, we take into account that the numbers of rulers and 
bandits decrease because of retirement and natural mor- 
tality. This loss is denoted by M(Y) for the bandits and 
R(Z) for the rulers. Taking all these factors into account 
one obtains differential equations of the form 

X = S(X) -C(X,Y) -T(X,Z), 

Y = rjC(X,Y)- L(Y,Z)~ M(Y), (1) 

Z = vC{X,Y) - R(Z) 

where r\ and v are constant factors. 



B. Normalization of the socio-economic model 

So far we have formulated a generalized version of the 
model proposed in [5|- In general the investigation of 
such a model would start with the computation of steady 
states. However, in our model we can not compute steady 
states with the chosen degree of generality. Instead, we 
assume that at least one steady state exists. This as- 
sumption is in general valid in a large parameter range. 
We denote the population densities in this steady state 
by A*, Y* and Z* respectively. 

Let us now consider the local asymptotic stability of 
the steady state under consideration. In a system of au- 
tonomous ordinary differential equations (ODEs), such 
as Eq. QJ, the stability of steady states depends on the 
eigenvalues of the corresponding Jacobian matrix • For 
any system ODEs of the form x = g(x) the Jacobian J 
can be computed as 



1,...,A. 



(2) 



In principle we could therefore compute the Jacobian of 
our system directly by applying this definition to the sys- 
tem from Eq. Q). However, the Jacobian obtained in 
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this way would depend explicitly on the unknown steady 
state. Because of this dependence the interpretation of 
any insights gained from such a Jacobian is difficult. In 
order to avoid this difficulty we study a normalized sys- 
tem. By means of normalization it possible to arrive at 
a system for which the Jacobian does not depend explic- 
itly on the unknown steady state. Instead, the unknown 
steady state only appears in certain parameters which 
are easily interpretable. 

In the following we use the abbreviated notation 



S* 
M* 



= S(X*),C* 
= T(X*,Z*),L* 
= M(Y*),R* 



= C(X*,Y*), 

= L(Y*,Z*), (3) 

= R{Z*). 



Using this notation, we define normalized state variables 

y = T^,z = — (4) 



X 

X* 



and normalized functions 

s(x) := 

c(x,y) := 

t(x, z) := 

l(y,z) := 

m(y) := 

r{z) := 



S{X*x) 

S* ' 
C(X*x,Y*y) 

C* 

T(X*x,Z*z) 

p* 

L(Y*y,Z*z) 
L* 

M(Y*y) 

M* ' 
R(Z*z) 



(5) 



R* 



Note that the normalized quantities have been defined in 
such a way that 

x* = y* = z* = s* = c* = t* = I* = m* = r* = 1 (6) 

holds. 

Let us now rewrite our model in terms of the normal- 
ized variables. We start by substituting the definitions 
Eq. J3J and Eq. © into Eq. (Q). We obtain 



S* C* T* 

T]C* L* M* 



(7) 



In order to simplify these equations we consider the 
system in the steady state. Because of Eq. © this yields 



S* 

X 7 ' 
rjC* 

Y* 
uC* 



C*_ 
~x* ' 

L* 
Y* 
R* 
~ TFT 



X*' 
M* 

Y* 



(8) 



Since Eq. JSJ contains only constants the equation has to 
hold even if the system is in a non-stationary state. It is 
therefore reasonable to define 



(9) 





S* 


c* 




a x = 


~X* 


~ ~x* 4 


- 


a y = 




L* 


M* 


Y* 


Y* 


' y* 


a z = 


vC* 


R* 




Z* 


- 





Note, that the third line of Eq. © allows us to replace 
all constant factors in the third line of Eq. Q by a z . In 
order to replace the other constant factors in the model 
in a similar way we define 



Ac := 
0y :-- 



1 C* 

cZx* 



1 T* 

^ - a y Y* 



(10) 



In the following we call a x , oty, ct z , 3 X and B y scale pa- 
rameters. The definitions Eq. Q and Eq. (jlt)|l effectively 
group the gain and loss terms together. As we will show 
below it is advantageous to define the scale parameters 
in this way since it results in easily interpretable param- 
eters. 

By applying Eq. © and Eq. ^| we can write Eq. (JJJ 

as 



x = a x (s(x) - B^c(x,y) - 3 x t(x,z)), 
V = oc y (c(x,y) -3yl(y,z) - B~ y m(y)), 
i = a z (c(x, y) - r(z)). 



(11) 



As a result of the normalization we have obtained a nor- 
malized model with the same structure as the original 
model. But, we know that the steady state under con- 
sideration is located at x* — y* = z* = 1 in the nor- 
malized model. One could argue that the normalization 
procedure is merely a transformation of difficulties: In 
order to perform the normalization we had to introduce 
the scale parameters that again depend on the unknown 
steady state. However, in the following we show that the 
scale parameters can be easily interpreted in the con- 
text of the model. That means, the values of the scale 
parameters can be determined by measurements or the- 
oretical reasoning rather than explicit computation. In 
this way the normalization procedure provides us with 
a set of parameters, depending on which the dynamics 
of the system can be studied. Let us now consider the 
interpretation of these parameters in more detail. 

From the way in which a x , a y and a z appear in 
Eq. (|llf> it can be guessed that these parameters de- 
note timescales. This guess is confirmed if one consid- 
ers Eq. @. In this equation a x is defined as the total 
production rate divided by the number of farmers. In 
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other words this parameter denotes the total per-capita 
production rate of the farmers in the steady state. The 
parameter a y is defined as the per-capita growth rate of 
the bandits. At the same time it has to be the per-capita 
rate at which bandits disappear (to prison, retirement or 
death). Consequently, we could write a y — l/r y , where 
r y is the typical length of a bandits career. In the same 
way a z is related to the typical length of a rulers career. 

The parameters /3 X and f3 y quantify the relative im- 
portance of certain processes. We have defined /3 X as the 
per-capita loss rate of farmers because of crime divided 
by the total per-capita loss rate. This means that /3 X is 
simply the fraction of the farmer's losses that occurs be- 
cause of crime. Likewise /3 X is the fraction of the farmer's 
losses that occurs because of taxes. In a similar way, the 
parameter (3 y denotes the fraction of the loss rate of ban- 
dits that occurs because of interaction with the rulers. In 
other words f3 y denotes the fraction of bandits that are 
eventually caught. At the same time [3 y is the fraction 
of bandits that escape the rulers till they retire or die 
because of natural mortality. 



C. The Jacobian of the socio-economic model 

So far we have derived a normalized model in which 
the location of the steady state under consideration is 
known. We have shown that the constant factors that 
appear in the normalization can be identified as mean- 
ingful parameters. Let us now compute the Jacobian of 
this normalized model. In the Jacobian all nonvanish- 
ing derivatives of the normalized functions with respect 
to the state variables appear. We consider these deriva- 
tives as additional parameters. In the following these pa- 
rameters are called exponent parameters. The exponent 
parameters are defined as 



malized function is simply 



9 ( \ 

-t(x,z) 
az 

9 ( \ 
az 



x=l 



*x 



x=l 



l v := 



x=l 



x=l 



(12) 



In order to understand the nature of the exponent pa- 
rameters it is useful to consider the effect of the normal- 
ization on a specific function. For instance in the specific 
model studied by 0, |f| the natural mortality rate of the 
bandits is modeled by the function 



M(Y) = AY, 



(13) 



where A is a constant parameter. According to our nor- 
malization procedure (see Eq. ©) the corresponding nor- 



m(y) 



AY 
~AY* 



(14) 



Therefore the exponent parameter is m y = 1 regardless 
of the value of A. In a similar way it can be shown that 
any specific function of the form 



M(Y) = AY q 



corresponds to an exponent parameter 



(15) 



(16) 



Since we have normalized all functions in the same way, 
this does not only hold for the mortality term M(Y) but 
for all functions in the model. For all processes that 
can be modeled with monomial expressions the expo- 
nent parameters are therefore identical to the exponent 
of the mononomial. If a process is not described by a 
mononomial then the corresponding exponent parame- 
ter still measures the nonlinearity of the process in the 
steady state. Consider for instance the specific function 



C(X,Y) 



AXY 

K + X 



(17) 



that is employed in [5| to describe crime. In this case our 
normalization procedure yields the exponent parameter 



1 



1 + X*/K 



(18) 



That means the parameter is almost one for X* <C K 
where C(X*, Y*) is almost linear in X*. But, it vanishes 
for X* >• K where C(X*,Y*) is almost constant in X*. 
For more examples of exponent parameters for specific 
functional forms see §. 

In this paper we will not assume that C(X, Y) or M(Y) 
have any specific functional form. But, we will use the in- 
sights gained from the consideration of specific funtional 
forms to interpret the exponent parameters. This in- 
terpretation reveals realistic ranges and values for these 
parameters. 

For instance, the parameter s x measures the nonlin- 
earity of the total production rate as a function of the 
number of farmers. In a country in which empty land of 
sufficient quality is available we can assume that the pro- 
ductivity increases linearly with the number of farmers. 
In this case we have s x ~ 1. However, if the production is 
not limited by the number of farmers, but by the lack of 
usable land then the production is a constant function of 
the number of farmers. This corresponds to s x = 0. We 
can therefore interpret s x as an indicator for the avail- 
ability of usable land in the steady state. 

Using the scale and exponent parameters the Jacobian 
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of the system can be written as 



£ X Px Px ^X Px Cy Px ^ Z \ 

C x Cy fty^y Py^^y fty^z I * 

C x Cy ~?"z / 

One of the parameters a x ,a y ,a z can always be set to 
one by means of timescale normalization. This leaves us 
with 13 parameters. For comparison the specific model 
proposed in p| contains 11 parameters. 



D. Computation of bifurcations 

The Jacobian which we have obtained in the previous 
section enables us to compute the stability of the steady 
states. In general the stability properties of dynami- 
cal systems change abruptly at critical parameter values 
which are known bifurcation points0, 0|. Let us focus 
on the bifurcations that are encountered as one parame- 
ter is changed. Of these local codimension-1 bifurcations 
only two types exist. The first type is characterised by 
the presence of a single zero eigenvalue of the Jacobian. 
In this paper we refer to these bifurcations collectively 
as bifurcations of saddle-node type. These bifurcations 
generally correspond to the emergence or destruction of 
equilibiria or the exchange of stability properties between 
two equilibria. Some examples of this type of bifurction 
are the saddle-node bifurcation, the transcritical bifur- 
cation or the pitchfork bifurcation. The bifurcations of 
saddle-node type can be found analytically by demand- 
ing that the determinant of the Jacobian vanishes. The 
other type of local codimension-1 bifurcation is the Hopf 
bifurcation. A Hopf bifurcation point is characterised 
by the presence of two purely imaginary eigenvalues of 
the Jacobian in the steady state. In this bifurcation the 
steady state becomes unstable as a limit cycle emerges or 
vanishes. In order to compute this bifurcation we use the 
method of resultants 01 ■ Although this method yields 
explicit analytical results the corresponding formulas are 
too lengthy to be presented here. Instead, the bifurca- 
tions are shown in a three-parameter bifurcation diagram 
in Fig. Q] 

In the figure the steady state under consideration is 
stable in the top-most volume of the parameter space. 
The stability of the steady state is lost in a Hopf bifurca- 
tion or in a bifurcation of saddle-node type. Since these 
bifurcations are of codimension one, the corresponding 
bifurcation points form hyper-surfaces in the parameter 
space. For the model under consideration the Hopf bifur- 
cation surface is of particular interest. As this bifurca- 
tion surface is crossed the steady state becomes unstable 
while a limit cycle emerges. This gives rise to the dynas- 
tic cycle that is observed in simulations. The fact that 
the bifurcation surface lies almost perpendicular to the 
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FIG. 1: Bifurcations in a society of farmers, bandits and 
rulers. Stable dynasties occur in the topmost volume of the 
parameter space shown. If the availability of usable land s x 
is increased the stability is lost as a Hopf bifurcation (dark 
surface) or a bifurcation of saddle-node type (light surface) 
is crossed. The Hopf bifurcation gives rise to oscillatory be- 
havior that is known as the dynastic cycle. The fraction of 
farmer's losses that occur because of crime (3 X and the frac- 
tion of the bandits that get eventually caught /3 y have only a 
weak effect on the stability of the dynasty. 



direction of s x indicates that this parameter is of greater 
importance for the stability of the system than /3 X or /3 y . 
Although such insights should be confirmed in a more 
careful mathematical investigation, they illustrate that 
bifurcation analysis of generalized models can be used to 
identify important parameters in a system. 

The scale and exponent parameters which we have de- 
fined above capture only the local behavior of the gen- 
eralized model. The bifurcation analysis which we can 
perform is therefore limited to the investigation of lo- 
cal bifurcations. Nevertheless, this analysis enables us to 
draw certain conclusions on the global dynamics of the 
generalized model. For example in Fig. ^ there is a line 
in which the Hopf bifurcation surface intersects the bifur- 
cation surface of saddle-node type. The points on such 
a line are codimension-2 Gavrilov-Guckenheimer bifur- 
cation points. The existence of these points implies that 
a region of complex (chaotic or quasiperiodic) dynamics 
exists generically in the model |9j- In fact, it is shown in 
[f| that chaotic dynamics appear in a specific version of 
the generalized model considered here. 

In summary we can say that the investigation of gen- 
eralized models allows us to generalize insights from spe- 
cific models. In particular, generalized models can be 
used to investigate the local dynamics (e.g. emergence 
of dynastic cycle from a Hopf bifurcation), to identify 
important parameters (e.g. availability of land s x ) and 
to draw certain conclusions on the global dynamics (e.g. 
the existence of complex dynamics). 
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III. A GENERALIZED MODEL OF COUPLED 
LASERS 

In the previous section we have shown that the local 
bifurcations of generalized models can be computed effi- 
ciently if the model is normalized in a certain way. This 
approach can be summarized in form of the following al- 
gorithm: 

1. Formulation of a generalized model 

2. Normalization 

• Define normalized variables and functions 

• Rewrite model in terms of the normalized 
quantities 

• Identify scale parameters 

3. Computation of the Jacobian 

• Define exponent parameters 

• Interpret the exponent parameters 

• Compute the entries of the Jacobian matrix 

4. Bifurcation analysis, . . . 

We argue that essentially the same procedure can be ap- 
plied to a large variety of different systems which can be 
described in the form of balance equations. In order to 
demonstrate this point we consider a system of coupled 
lasers. Although the scientific background is entirely dif- 
ferent, it turns out that the laser model can be treated in 
the same way as the socio-economic model. For the anal- 
ysis presented here the most striking difference between 
the physical and the socio-economic model is the amount 
of available information. In the socio-economic model 
all processes are difficult to quantify and are therefore 
candidates for generalization. By contrast, the underly- 
ing quantum mechanical foundations of laser physics are 
well understood. Many of the processes in the model can 
therefore be quantified with a high degree of certainty. 
Nevertheless, it can be useful to generalize certain pro- 
cesses which can be in principle quantified, but depend 
on many details of the specific system under considera- 
tion. In the following we show that the algorithm out- 
lined above can be applied in such a partially generalized 
system. 



A. Formulation of a model for coupled lasers 

Models from laser physics are among the classical ap- 
plications of the theory of dynamical systems fTTL fl~3 . fl3| . 
In particular coupled lasers are known to exhibit inter- 
esting dynamics. The investigation presented here is in- 
spired by the theoretical and experimental results re- 
ported in [I I Il5 | . In this work it is shown that the 



dynamics of two coupled single-mode lasers can be de- 
scribed by equations of the form 



Ai = eiAid + rid + k^A 1 A 2 cos(<5) - L 1 (A 1 ) 
A 2 = e 2 A 2 G 2 + r 2 G 2 + kV 'AiA 2 cos(tf) - L 2 (A 2 ) 
G\ = pi — siGi — r\G\ — eiA\G\ 
G 2 = p 2 — s 2 G 2 — r 2 G 2 ~ e 2 A 2 G 2 

In order to discuss these equations, let us use the index 
i = {1,2} to refer to laser 1 and 2 respectively. In the 
equations the intensity of laser i is described by the vari- 
able Ai while the gain of each laser is described by the 
variable Gi. The variable 8 denotes the phase difference 
between laser 1 and laser 2. 

The equations that governs the time evolution of Gi, 
the gain of laser i, contains four terms. The first term pi 
represents a constant pumping of the laser. The second 
term SiGi describes the decrease of the gain because of 
spontaneous emission of photons which do not go into 
the laser mode. The term r^G^ represents spontaneous 
emission into the laser mode. The final term eiAiGi cor- 
responds to stimulated emission. The quantities Pi, s,-, 
Ti and &i that appear in these equations are assumed to 
be constants. 

Spontaneous or stimulated emission of photons into the 
laser mode of laser i increases the intensity Aj. This is de- 
scribed by the first two terms in the intensity equations. 
The third terms in these equations models the optical 
coupling of the lasers. The specific functional form of 
these terms as well as the equation that governs the evo- 
lution of the phase difference are results of Stratonovich 
calculus 01 • In these equations the constant n describes 
the strength of the coupling. The parameter r c is the cav- 
ity round trip time and i5o is the difference of the phase 
velocity of the lasers without coupling. 

The only general function that appears in the model 
equations is Li(Ai). This function is used to describe 
the loss of laser intensity. Photons can be lost from the 
laser mode because of the absorbtion in the medium or 
because of escape from the resonator. Nonlinear optical 
effects like multi-photon absorbtion introduce nonlinear- 
ities in the intensity dependence of absorbtion inside the 
resonator [l6|. Other nonlinear effects, like for instance 
temperature dependence or intensity dependence of the 
refractive index can effect the goodness of the resonator. 
As a consequence the losses of intensity can depend non- 
linearly on the intensity. While all of these effects can be 
quantified for a specific laser medium, this would yield a 
complicated model that describes only one very specific 
type of laser. In many models studied in literature it is 
instead assumed that the losses depend linearly on the 
intensity of the laser. 



7 



B. Normalization of the physical model 



identify the scale parameters 



Let us now normalize the generalized model formu- 
lated above. We start by assuming that there is at least 
one steady state. The dynamical variables in the steady 
state are denoted by A\* , A 2 * , Gi*, G 2 * and S* respec- 
tively. Furthermore, we denote the values of the gen- 
eralized functions in the steady state by L*. As in the 
previous example we define the normalized state variables 



A>* 
Gj_ 

G/ 



(21) 



and the normalized function 



k(a.i) := 



(22) 



Normalization of the variable S is not necessary. We 
could argue that S is already an easily observable and 
interpretable variable, therefore the steady state value 
5* can be treated as a scale parameter. However in the 
following it will turn out that this reasoning is not nec- 
essary. 

In the next step we rewrite the model in terms of the 
normalized variables. By substitution of Eq. 121|) and 
Eq. into Eq. we obtain 



eiAi*Gi* , rid* . k\J A x * a\A 2 a 2 



A 



A, 



A, 



cos(<5) 



Pi SiGi TiGi _ &iAi Gi 
Gi Gi 



tt9i n * yi 



Gi 



9i 



Gi 



r c \ V A 2 a 2 y A 1 ai J 

(23) 

where we have used the abbrieviating notation L* — 
Li(A*), In order to identify suitable scale parameters 
we consider the system in the steady state. This yields 



= ad + — — + cos((5 ) - — 

Ai J±i J±i 

= 77T - Si - n- e l A i 



(24) 

In the socio-economic model studied in the previous sec- 
tion we have seen that it is advantageous to group gain 
and loss terms together. In this way we obtain scale 
parameters that correspond to interpretable time scales. 
Grouping the terms in the laser model in this way we can 



a, 



Pi ■= 



eiGi 

Pi 
G,* 



r i Gi 

~a7 



k v / A 1 *A 2 



A 4 



■ cos((T) + -hr 



i.Ai 



. (25) 

Note that we have treated the coupling term in the inten- 
sity equation as a loss term. This decision is motivated 
by the fact that frequency locking in experiments is gen- 
erally observed to be anti-phase [l5j . In this case 5* s» 7r 
and cos(5*) < 0. The analysis presented below will show 
that it is generally reasonable to consider the coupling 
term as a loss term. Our mathematical treatment of the 
system does not depend critically on the negativity of the 
coupling term. However, if the coupling term were posi- 
tive it would be advisable to define the scale parameters 
in a different way for the sake of interpretation. 

The way in which the scale parameters have been de- 
fined above allows us to identify the parameter as 
the characteristic timescale of the intensity equation if 
7r/2 < S* < 3ir/2. In this case is the per-capita 
(per-photon) loss rate of intensity of photons in the laser 
mode. In other words this ai is the inverse of the typical 
lifetime of a photon in the laser mode of laser i. Likewise 
Pi is the characteristic timescale of the gain equation or 
the inverse of the typical lifetime of atoms in the upper 
laser level. 

In order to eliminate the constant factors from Eq. 123|) 
we need to define additional scale variables which mea- 
sure the contribution of the individual terms to the total 
gain and loss rates. We define the additional scale pa- 
rameters 

- l^ll 

Xi ■ a * J 

a, Ai 

a i := ir( r i + e i A i*), 
1 

It ■= —eiGi* 
and the complementary parameters 



(26) 



, 1 M 

Xi = 1-Xi = cos{d ), 

a t Ai 



1 

o-i = 1 - &i = -5- Si, 

1 TiGi* 

7i = 1 - 7» = 



(27) 



ai Ai 



Note that 



Finally we define 



PiO; 



■dAi* ,74 = 



Pidi 



A\ 

P : = US- 



(28) 



(29) 



We can interpret the scale parameters as follows. The pa- 
rameter Xi denotes the fraction of the losses of intensity 



8 



in laser i that occur because of absorbtion or losses from 
the resonator. The complementary parameter Xi is the 
fraction of the losses of laser i that occur because of the 
coupling to the other laser. The parameter describes 
the fraction of the total loss of gain of laser i that corre- 
sponds to the production of photons in the laser mode. In 
a similar way the parameter ji denotes the contribution 
of stimulated emission to this fraction. In other words, 
7, denotes the fraction of photons in the laser mode of 
laser i that were created by spontaneous emission. 

The scale parameters allow us to write the model equa- 
tions as 



, - - / cos(<5) 

CMTiOiffi + 7*0, - XiV a i a 2 TTTT 

COS((T) 

A (1 - ~ Vi(cim + 7iOi5i)) 



5 = 




sin(<5) — Sq 



(30) 



C. Computation of the Jacobian 



Having completed the normalization we can now com- 
pute the Jacobian of the model. We start by defining the 
exponent parameter 



-k{cti) 



(31) 



Qi = l 



Using this parameter we can write the Jacobian in the 
steady state as 



7i 



a 2 



02 



k/t c / 



V 0.5(p 



1/2 



0-5xi - XiMi 
-0.5X2 
— ffi7i 


_ n-V2 



72 



-0.5xi 

0.5X2 - X2M2 





)sin(<T) 0.5(p 



-1/2 



-0-272 
-P 1/2 ) 



sin(5* 



(p^ + p-^cosiS*) J 



Xi taa(<S*) 
X2tan(5*) 





l/2\ 



(32) 



Let us assume that both lasers operate with same inten- 
sity in the steady state. This implies p = 1. In this case 
all elements of the fifth row of the Jacobian vanish ex- 
cept for Js : 5. Therefore we can immediately identify the 
eigenvalue A5 = = 2cos(<5*). The existence of this 
eigenvalue proves that stable steady states are only pos- 
sible for 7r/2 < S* < 3ir/2. In this interval the stability of 
steady states is determined by the other four eigenvalues 
and therefore by the reduced Jacobian 



7i 







/ «i 


\ 








Jrcd 




V 


a 2 

Pi 

P2 I 


X 






0.5xi - 




-0.5xi 




1 





-0.5X2 




72 


- 0.5x2 - X2M2 





1 


-ci7i 











— I 












— 0272 







-1 



(33) 



D. Computation of Bifurcations 

Let us now study the bifurcations of the generalized 
laser model based on the reduced Jacobian. We simplify 
the reduced Jacobian further by assuming that the two 



lasers in the system are described by identical parameters 
(ax = a 2 =: a, (3 X = /3 2 =■ P, 7i = 72 =: 1, Xi = X2 =: 
X, fJ-i — fJ-2 —'■ fJ- and (Ti = (7 2 =: a). This assumption is 
valid for many coupled laser systems that are studied in 
experiments. Furthermore we set a = 1 by means of time 
normalization. As in the previous example we apply the 
method of resultants to compute the local bifurcations 
[Tof . This reveals the bifurcation hypersurfaces of saddle- 
node type 



SI : /i = 



7(0- -1) 



for a < i 



7(<t-1) + 1-X 



S2 : fi = — — - for a < £ 

X 1 

and the Hopf bifurcation hyper-surfaces 

7 — 6 



(34) 



HI 



(J- 



H2 : fj, 



6 + X-l 



X 



for a < S- 

1 

for a < - 

7 



(35) 



These bifurcation surfaces are shown in Fig. The 
steady state under consideration is stable in the top-most 
volume in the front of the left diagram. The right dia- 
gram shows the bifurcation surfaces from a different per- 
spective so that the stable area is now in the back of the 
diagram. The stability is lost if one of the Hopf or saddle- 
node bifurcation surfaces is crossed. Note, that for a < 
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FIG. 2: (Color) Two different perspectives on a three pa- 
rameter bifurcation diagram of the coupled laser system. 
The steady state under consideration is stable in the vol- 
ume that is located at high values of \ and fi. The di- 
agrams show the two saddle-node bifurcation surfaces SI 
(dark blue) and S2 (light blue) and the two Hopf bifurca- 
tion surfafaces HI (red) and H2 (green). At the intersection 
of these surfaces codimension-2 Takens-Bogdenov, Gavrilov- 
Guckenheimer and double Hopf bifurcations are formed. All 
codimension-2 lines meet in a codimension-4 quadruple point 
bifurcation. While the codimension-1 indicate the borders of 
the region of local stability, the bifurcations of higher codi- 
mension enable us to draw certain conclusions on the global 
dynamics. 



6/7 both Hopf bifurcation surfaces are independent of a. 
At er = 6/7 both surfaces end simultaneously in Takens- 
Bogdanov bifurcation lines. The Takens-Bogdanov line 
of the Hopf bifurcation surface HI is formed as it meets 
the bifurcation surface SI. At the same value of a the 
surface H2 ends in a Takens-Bogdanov bifurcation line 
as it meets S2. In addition to the Takens-Bogdanov lines 
there are two other codimcnsion-2 bifurcation lines. A 
Gavrilov-Guckenheimer bifurcation is formed as H2 in- 
tersects SI. Finally, HI and H2 meet in a double Hopf 
bifurcation line at \ — 1. All of these codimension-2 bi- 
furcation lines meet in a single point at \ — 1, a = 6/7. 
This point is a bifurcation point of codimcnsion-4. 

As in the previous example the existence of bifurca- 



tions of higher codimension enables us to draw conclu- 
sions on the global dynamics of the system. The pres- 
ence of Takens-Bogdanov bifurcations proves that a ho- 
moclinic bifurcation has to exist. This bifurcation corre- 
sponds to the formation of a homoclinic loop that con- 
nects one saddle to itself. Homoclinic loops are known 
to cause bursting behavior. In laser experiments such 
dynamics were observed in |l5l IT?! ] . Furthermore homo- 
clinic bifurcations play an important role in the forma- 
tion of Shilnikov chaos [l^. In the previous section we 
have already mentioned that the Gavrilov-Guckenheimer 
bifurcation indicates the existence of a region of complex 
(at least quasiperiodic) dynamics. In a similar way the 
double Hopf bifurcation indicates the generic existence 
of a region in which chaos can be, at least transiently, 
observed 0. 

The codimension-2 bifurcation lines do not only show 
that certain global dynamics do generically exist in the 
system. They also point to the parameter region where 
such dynamics are likely. For instance chaotic dynamics 
is likely close to the double Hopf bifurcation line. In this 
way the relatively simple computation of codimension-2 
bifurcations can be used to locate interesting regions in 
parameter space. In these regions more time consum- 
ing techniques, like for instance the computation of Lya- 
punov exponents in specific models can then be applied. 



IV. A GENERAL FOOD WEB MODEL 

In the previous sections we have applied the proposed 
approach to a generalized social model and to a partially 
generalized model of coupled lasers. In this section the 
same procedure is applied to a generalized model, which 
describes the dynamics of ecological food webs. The long 
term behavior and in particular the stability of food webs 
is investigated in a l arge number of theoretical and obser- 
vational studies [h1I201|21| . In the following we compute 
the Jacobian of a generalized food web model. Since this 
model is more involved than the ones considered in the 
previous sections some additional difficulties arise. Nev- 
ertheless, the Jacobian can be constructed by following 
the algorithm outlined in Sec. IIIII 



A. Formulation of the ecological model 

Let us consider a food web of N ecological populations. 
We denote the size of these populations with dynamical 
variables X±, . . . , Xn . For the sake of generality we do 
not distinguish between producer and consumer popula- 
tions. Every population in the model can grow either by 
primary production (e.g. photosynthesis) or by consump- 
tion of prey (predation). The size of a given population 
can decrease because of predation by others or because of 
other causes of mortality (e.g. natural aging, diseases). 
We denote the primary production of population n by the 
function S n (X n ). The growth of population n by preda- 



10 



tion on others is described by the functional response 
F n (Xi, . . . , Xjy). The loss inflicted by a higher predator 
m is described by the function L mtn {X\, . . . , Xn). The 
mortality of species n is denoted as M n (X n ). Finally, we 
introduce the parameter rj which denotes the portion of 
consumed prey biomass that is converted into predator 
biomass. This yields the model 

X n = S n (X n )+r ]n F n (X 1 ,...,X N )-M n (X n ) 

A (36) 



L m , n (Xi, . . . , Xn) 



where n = 1, . . . , N. 

In order to describe natural food webs in a realistic 
way we have to take an additional constraint into ac- 
count. The loss rate of population n because of preda- 
tion by population m is not unrelated to the growth rate 
F m of population m. In order to quantify this relation- 
ship we introduce the auxilliary variable T m (X\, . . . , Xn) 
which denotes the total amount of prey that is available 
to species m. It is reasonable to write T m {X\, . . . ,Xn) 
as 



/v 



T m {X\, . . . ,X[f) — C rn ,n{X n ), 



(37) 



where C m ,n{X n ) are auxilliary variables that describe the 
contribution of the population n to the total amount of 
food available to population m. The value of C m . n {X n ) 
depends not only on the density of population n but also 
the preference of the predators for specific prey and on 
the success rate of attacks [j^ • 

We can now write the amount of prey consumed by 
population m as 



F m (Xi, . . . , X N ) — F m (T m , X m ). 



(38) 



Since species n contributes a fraction 
C m ,n{X n )/T m (Xi, . . . ,Xn) to the total amount of 
prey available to species m it is reasonable to assume 
that it contributes the same fraction to the consumption 
rate. That means, 



L m ,n{Xi, . . . , Xn) 



Cm : n(X n ) 

T m (X\, . . . , Xn) 



Fm {F m , X m ) 



(39) 



B. Normalization of the ecological model 



So far we have formulated a model of a TV-species food 
web. In contrast to the models considered in the previous 
section we had to impose additional constaints in form of 
the algebraic equation Eq. I|39|) . In the following we apply 
our normalization procedure to normalize these algebraic 
equations in addition to the ODEs. 

Like in the previous sections we use asterics to indicate 
the values of functions and variables in the steady state. 



Using this notation, we define normalized state variables 

*n = ^ (40) 

normalized auxilliary variables 

_ T m (Xi*x 1 , . . . ,X N *x N ) 

• rp * 3 

± 171 

Cm , n (-^n %n) 

Cm,? i ■ — * 
(- / m,n 

and normalized functions 



(41) 



s n {x n ') 
f n (xi,...,x N ) := 

lm,n( x li • ■ • j x n) '•- 



M n (X n *x n ) 
F n (Xi*xi, . . . , X N *x N ) 

L m ,n{Xi*Xi, . . . , Xn*Xn) 

f * ' 



(42) 



As the second step of the normalization we rewrite the 
model in terms of the normalized quantities. We start by 
substituting the definitions Eqs. l|10 ]) -l|i2 )l into Eq. 
This yields 

Xn "TT * S n \X n J -\- — j J n yX\ , . . . , Xn ) 

M* N T * 

T m n {X n ) - > - ; l m ,n(Xl, ■ ■ ■ , Xn) 

— \ -&n 

m—1 



x„ 



(43) 

In the next step we have to define suitable scale variables. 
We start with the time scales 



Sn* , VnFn* _ M* 



N 



T * 



V V V / j y 

A-n -A-n ***-n -, -™-n 

m—1 



(44) 



These scale parameters quantify the rate of biomass flow 
in the steady state. The relative contributions of the 
different processes to the biomass flow can be described 
by 



Pn := 



1 y n F n 

On X n * 



Pn ■= 1 - Pn = 



N 



1 s n * 

X* 1 
n 



1 _ _ T * 



a n ' X 

m—1 



(45) 



a n := 1 - a n = 



1 L r . 



1 M* 

Ot n X n * 



Before we discuss the interpretation of these parame- 
ters in more detail let us examine the algebraic equations 
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in the model. We start with Eq. i|39[l . In this case the 
substitution of the normalized quantities yields 



lm,n(Xli ■ ■ ■ i %N) 



C *F * r 

rp * T * i 

Cm ' n f (f r ) 



fm (fm i X ra ) 



(46) 



where we have used Eq. (|39|) in the second step. 

The last equation which we have to consider is Eq. (p?7j) . 
In this case we obtain 



N 



km ^ ^ rr- ''ni.i!- 



n=l 



Tn 



(47) 



In contrast to Eq. (|46|) the constant factors do not cancel 
out in this equation. It is therefore useful to define the 
additional scale parameters 



Xr. 



(48) 



We can now write the food web model as 

&n ^■ni^Pn^n^n) ~\~ Pnfn\^n^n) ^"n^nO^n) f a(\\ 
_ R 1 I W ^ i 



for n — 1 . . . N and 



^m,n(^l: • • • ; *£./v) 



" fm (^m : •Em) 1 (^^) 



the growth terms, a n and a n denote the relative impor- 
tance of loss terms. The parameter a n is defined as the 
fraction of the total loss rate of population n that occurs 
because of predation by explicitly modelled predators. 
The fraction of losses that is caused by all other sources 
of mortality (e.g. natural mortality, diseases, predation 
by predators that are not explicitly modeled) is given by 
a n . We can therefore say, that the value of a n indicates 
the predation pressure on species n. 

The parameter (3 m ^ n is defined as the per-capita loss 
rate of population n because of predation by population 
m divided by the total loss rate of population n because 
of predation a n a n . This means that /? mj „ denotes the 
relative weight of the contribution of species m to preda- 
tive loss rate of species n. If, for example, population 
1 is the only population that preys upon population 2 
then the corresponding parameter is /3i : 2 = 1- However, 
if population 1 causes only 50% of the losses while the 
other half is caused by predation by cannibalistic preda- 
tion of population 2 upon itself then we find = 0.5, 
#2,2 = 0.5. 

Finally, the parameter Xm,n measures the relative con- 
tribution of population n to the total amount of food that 
is available to species m. Note that, Xn.m (3 n ,m- For 
example predation by species 1 may only contribute, say 
l/10th of the total predative loss of species 2. But, at the 
same time species 1 may be the only species which is con- 
sumed by species 2. In this situation we have (3\^ = 0.1 
but xi,2 = 1- 



(51) 

n=l 

As a result of the normalization we have identified a set 
of scale parameters. Let us now discuss these parameters 
in more detail. 

Based on our experience from the previous section we 
can identify the parameters a,\, . . . , czn immediately as 
the characteristic timescales of the system. The param- 
eter a n denotes the per-capita birth and death rate of 
population n in the steady state. In other words, a n 
is the inverse of the typical lifetime of individuals from 
population n. 

The parameter p n is defined as the per-capita biomass 
gain of population n divided by the total per-capita 
growth rate a n . In other words, this parameter de- 
scribes what fraction of the growth rate of population 
n is gained by predation. The complementary parame- 
ter p n denotes the fraction of the growth that originates 
from primary production. If population n consists of 
primary producers the corresponding parameter values 
are p n = 0, p n = 1. A population of consumers has 
p n = l,p n = 0. Some organisms are known which are 
capable of primary production and predation, for these 
p n and p n can have fractional values between zero and 
one. 

The parameters a n and a n are very similar to p n and 
p n - While p n and p n describe the relativ importance of 



C. Computation of the Jacobian 

Like in the previous examples we start the computation 
of the Jacobian by defining a set of exponent parameters 
that describe the nonlinearity of the ecological processes 
in the steady state 
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X n 



(52) 



To a large extend these definitions follow the same 
scheme that we have applied in the previous sections. 
Note however that t n * appears in the definition of tpn- 
In order to understand the reason for this definition we 
have to be aware of the fact that cannibalism can occur. 
Consequently, the amount of food t n that is available to 
the predator depend on x n itself. The definition 

given above ensures that tp n denotes only the nonlinearity 
of the predator-dependence of the response function. 

The derivation of the Jacobian for the general food 
web model is analogous to the previous examples. For 
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this reason the calculation is presented in appendix D As 
a result we find that the non-diagonal elements of the 
Jacobian can be written as 

Jn,i — &n(Pn^fnXn,i^n,i 0~7l\fii,nlfti (53) 

(7m - l)Xm,«)) 

and the diagonal elements as 

Ji,i = Oti{Pi<t>i + Pi(liXi,i^i,i + 

-aifii- <Ji(Pi,iipi+ (54) 

~t~ X-^m— 1 /^m,i'^m,i((0 / m l)Xm,i "I - 1))) 

The parameter </>„ that appears in these equations de- 
scribes the nonlincarity of the primary production of 
species n. In an environment where nutrient are abun- 
dant and no other limiting factors exist, it is reasonable to 
assume that the primary production S n is proportional to 
the number of primary producers X n . This corresponds 
to a n = 1. On the other hand, if nutrients are scarce then 
the primary production is limited by the nutrient supply. 
In this case the primary production is a constant func- 
tion of the number of primary producers and a n — 0. In 
a realistic situation o~ n has a value between and 1. This 
value can be interpreted as an indicator of the nutrient 
availability in the system. 

The parameter fj, n denotes the nonlinearity of the mor- 
tality rate of population n. For top predators this param- 
eter is often called exponent of closure in the ecological 
literature. More generally speaking, we call fi n exponent 
of mortality. In most food web and food chain models 
this parameter is assumed to be either one or two. How- 
ever, it has been shown that fractional values between 
one and two may describe natural systems in a more re- 
alistic way. This is discussed in detail in [23|. 

The parameter j n denotes the nonlinearity of the pre- 
dation rate of population n with respect to the prey den- 
sity. If the prey is abundant then the predation rate is 
not limited by prey availability and j n — 0. However, 
in most natural systems prey is relatively scarce. In this 
case the value of 7„ can depend on the feeding strategy 
of population n. In the limit of scarce prey the value of 
7„ approaches one for predators of Holling type II and up 
to two for predators of Holling type III. This is discussed 
in more detail in 0. 

The parameter ip n describes the dependence of the pre- 
dation rate of population n on the population density of n 
itself. In most models it is assumed that predators of the 
same species feed independently without direct interfer- 
ence. In this case the predation rate is a linear function of 
the predator density and ip n = 1. However, intraspecific 
competition or social interactions can introduce nonlin- 
earities in the predator dependence of the predation rate. 
If these effects are taken into account we find < ip n < 1. 
An extreme example are the so-called ratio dependent 
models [24j . These models correspond to %j) n = 1 — j n . 

Finally, the parameter A m .„ describes the nonlinear- 
ity of the contribution of population m to the diet of 
population n. In the simplest case we can assume that 



the predator does not distuinguish actively between dif- 
ferent prey types. An example of such a predator is an 
aquatic filtration feeder. In this case the contribution 
of a suitable prey population m is proportional to the 
size of the population m and X m ,n = 1- However, other 
predators attack their prey individually. These predators 
may adapt their strategies to a specific prey population if 
they have sufficient practice. In this case the success rate 
of attacks as well as the enounter rate increase linearly 
with the density of the prey population. In effect that 
leads to a quadratic dependence and therefore A m ,„ = 2. 
Another type of behavior is exhibited by predators that 
need to choose their prey actively in order to obtain all 
nutrients that are necessary for growth. As an extreme 
example we could imagine a predator that requires ex- 
actly the right composition of its diet. In this case we 
find X m ,n = for all prey populations m but the limiting 
one. 



D. Stability of generalized food webs 

In Sec. Inland Sec. lIIII we have used bifurcation analysis 
to study the Jacobians of the normalized models. For the 
investigation of food webs this analysis is also possible. 
In fact, bifurcation analysis has been used in 0,0, El to 
study the dynamics of generalized food chains. The gen- 
eralized food chains considered in these papers are par- 
ticular examples of the more general class of food webs 
considered here. The bifurcation analysis of other exam- 
ples of generalized food webs can reveal new ecological 
insights. For instance we can consider bifurcation dia- 
grams that correspond to different food web geometries. 
In this way one can investigate which food webs are char- 
acterized by a high local stability. Furthermore we can 
investigate the impact of certain ecological mechanisms 
on the stability or use higher codimension bifurcations to 
draw conclusions on the existence of complex dynamics. 
However, from the methodological point of view, these 
investigations are to a large extent parallel to the previ- 
ous examples. Let us therefore describe a different way 
in which the derived Jacobian can be applied. 

In order to investigate the stability of ecological sys- 
tems random matrix models arc frequently considered 
[25| . In this context the Jacobian of a food web is mod- 
eled by a random matrix. The eigenvalues of a large set 
of these randomly created matrices are then computed. 
In this way one can for instance obtain a relationship 
between the stability of a random food web and the con- 
nectivity (often measured in terms of the number of non- 
diagonal elements of the matrix). While these analyses 
have provided many interesting results it has often been 
argued that not all of the randomly created matrices are 
biologically reasonable |26| . 

The Jacobian of a generalized food web of given size 
can be studied in a similar way. By choosing the scale 
and exponent parameters randomly from suitable distri- 
butions one can create a set of random but realistic Ja- 
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cobians. These Jacobians can then be studied by the 
numerical computation of eigenvalues. The advantage of 
this approach is that random matrices with certain eco- 
logical properties can be created. In this way one can 
for instance investigate the impact of ecological effects 
like intraspecific competition or cannibalism on the sta- 
bility of the system. Although these investigations are 
certainly promising they require a much more detailed 
analysis and are therefore beyond the scope of the cur- 
rent paper. 



V. DISCUSSION 

In this paper we propose an approach to the inves- 
tigation of generalized models which can be formulated 
in form of balance equations. In these equations gain 
and loss terms determine the time evolution of the state 
variables of the system under consideration. The pro- 
posed normalization procedure enables the researcher to 
consider a system from an abstract point of view. This 
abstractness is caused by the way in which variables and 
functions are normalized. The normalization procedure 
separates the necessary information about the steady 
state (the scale parameters) from the information about 
the nonlinearities in the model (the exponent parame- 
ters). In this way it becomes possible to study the effect 
of the nonlinearities independently from the location of 
the steady state. Apart from enabling us to investigate 
generalized models, this separation of scale and exponent 
variables has other advantages. For instance the elements 
of the Jacobians that are derived in this way are gener- 
ally relatively simple. More importantly, it is in almost 
all models possible to identify the parameters in such a 
way that they are easily interpretable in the context of 
the model. In many cases the general parameters are 
more directly accessible to measurements than the spe- 
cific parameters that are normally used in modeling. In 
it has been shown that the identification of the general 
parameters can in itself provide new insights. 

The basic idea behind the proposed approach is to con- 
sider the bifurcations as functions of a new set of param- 
eters which only encode as much information as is needed 
to construct the Jacobian. In this way we avoid the com- 
putation of steady states, which is in many cases more 
complicated than the computation of bifurcations once 
the Jacobian is known. This enables us to compute the 
local bifurcation surfaces analytically (using computer al- 
gebra systems) as a function of the scale and exponent 
parameters. As a result one is able to investigate the 
role of the nonlinearities in the gain and loss terms inde- 
pendently from the magnitude of these terms. In other 
words, one can study, with a high degree of generality, 
which functional forms of the gain and loss terms yields 
which dynamical behavior. Thus, this approach provides 
a bridge between modeling and stability analysis in dy- 
namical systems theory. 

The focus of the bifurcation analysis is the compu- 



tation of local codimension-1 bifurcation of the normal- 
ized steady state. But, as a byproduct of this analysis 
we obtain the local bifurcations of higher codimension, 
which can reveal many insights about the dynamics in 
the neighborhood of these bifurcations. For instance the 
presence of a codimension-2 double Hopf bifurcation is 
an indicator for chaotic dynamics. Although our analy- 
sis is purely local, it can therefore reveal insights on the 
global dynamics of the system. In particular it can indi- 
cate regions in parameter space where complex dynamics 
can be expected. This insight can be used to prove the 
generic existence of chaotic parameter regions in general- 
ized models or to locate interesting regions for subsequent 
numerical investigation. 

At present there are only few results on the implica- 
tions of local bifurcations of codimension larger than two. 
However, the opinion has often been expressed that this 
line of mathematical research would benefit from more 
examples from physical applications pj. The approach 
proposed here provides physicist with a tool for finding 
these examples. In return, new mathematical results on 
the implications of bifurcations of higher codimension 
can be applied directly to the investigation of general- 
ized models of physical systems. In this way, the inves- 
tigation of generalized models may lead to more efficient 
transport of mathematical discoveries into physical ap- 
plications and back. 

Our approach to the investigation of generalized mod- 
els is purely determistic. Noise and stochastic environ- 
mental fluctuations can not be accomodated adequately 
in the models considered here. However, the aim of the 
proposed approach is not to predict the time evolution of 
a given system, but rather to gain a general understand- 
ing of the functioning of a class of models. For systems 
with a low or moderate level of noise, the investigation of 
deterministic generalized models can still be useful, since 
it is unlikely that the presense of noise alters the local 
bifurcation structure qualitatively in this case. 

The abstract nature of the proposed analysis becomes 
apparent if one considers that we did not have to specify 
which steady state we consider. In a system in which 
multiple steady states exist the derived bifurcation dia- 
gram does therefore describe all of these steady states 
simulateously. This insight can for instance be used to 
prove general statements of the type: "There can not be 
a stable steady state with . . . " . 

In this paper, the main tool for the investigation of the 
derived Jacobians was bifurcation analysis. However, we 
have already remarked in our final example that the Ja- 
cobians can also be used in other ways. In particular 
large systems may contain too many parameters to be 
analysed efficiently with bifurcation diagrams. For these 
systems our method can be used in the random-matrix- 
like context, which is described in Sec. II VI In this light, 
the proposed method can be seen as a link between the 
very general random matrix approach and specific mod- 
els. 

Finally, let us remark that the approach presented here 
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is not restricted to systems of ordinary differential equa- 
tions. The algorithm which we have outlined in Sec. 11111 
can be applied in the same way to compute the bifurca- 
tions of discrete time maps. Moreover, it can be used in 
certain systems of partial differential equations to com- 
pute local bifurcations of homogenous solutions including 
Turing bifurcations [27j ■ 
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Using this result we compute 



^m,n (^1 ) • • * ) %N ) 



dx 



where 




Xm.i^m.i 7mXm,z^m,'i ^ 
(7m l)Xfi,i^m,i ~t~ Ql- 

(A.4) 



for n =/= i,m =/= i 
for n = i,m =/= i 
for n =/= i,m = i 
for n = i.m = i 



(A.5) 



Let us now compute the derivatives of the sum in 
Eq. (gSJ) . We obtain 



APPENDIX: COMPUTATION OF THE 
JACOBIAN OF THE GENERAL FOOD WEB 

In this appendix we present the explicit comptation 
of the Jacobian of the food web model. We start by 
considering 



Q x . ^2trn— 1 Prn,nlm,n{%li ■ ■ ■ > ^iv)) 



■7^T-t m (xi, . . . , X]\[) 
= "dx~ £n=1 Xm.nCm.n(Xn) 

x— x* 

We write the derivative of the predation rate as 



(A.I) 



+ 



(A.2) 



This yields 



d_ 

dx, 



fn {pn i 



lnXn,i^n,i for 1 ^ 71 

TnXn,iK,i + ipi for i = n 



(A.3) 



where 



for i = n 
for i 7^ n 



(A.6) 



(A.7) 



Using this result we can now write the non-diagonal ele- 
ments as 



n I, . i 1 



and the diagonal elements as 

Ji,i = Ctiip^i + Pi{"tiXi,i\,i + i>i) 

-ZiPi - (TiiPijipi (A. 9) 

~t~ ^^m—l $rn,i ^m,i ( (7m l)Xm,i ~t~ 1))) 
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